/**
 * @file sfmt.h
 *
 * SFMT random number generator
 *
 * The code is adapted from dSFMT Library by Mutsuo Saito and Makoto Matsumoto
 * (see http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/)
 *
 * @author Dahua Lin
 */

#ifdef _MSC_VER
#pragma once
#endif

#ifndef GMLPP_SFMT_H_
#define GMLPP_SFMT_H_

#include <gmlpp/base/basic_defs.h>
#include <gmlpp/base/simd.h>
#include <ctime>


namespace gml
{

	struct SFMTParamsCommon
	{
		static const uint64_t LOW_MASK = 0x000FFFFFFFFFFFFFULL;
		static const uint64_t HIGH_CONST = 0x3FF0000000000000ULL;
		static const int SR = 12;
	};

	template<uint32_t MExp> struct SFMTParams;

	template<>
	struct SFMTParams<19937> : public SFMTParamsCommon
	{
		static const uint32_t MEXP = 19937;

		static const uint32_t N = ((MEXP - 128) / 104 + 1);
		static const uint32_t N32 = N * 4;
		static const uint32_t N64 = N * 2;

		static const int POS1 = 117;
		static const int SL1 = 19;
		static const uint64_t MSK1 = 0x000ffafffffffb3fULL;
		static const uint64_t MSK2 = 0x000ffdfffc90fffdULL;
		static const uint32_t MSK32_1 = 0x000ffaffU;
		static const uint32_t MSK32_2 = 0xfffffb3fU;
		static const uint32_t MSK32_3 = 0x000ffdffU;
		static const uint32_t MSK32_4 = 0xfc90fffdU;
		static const uint64_t FIX1 = 0x90014964b32f4329ULL;
		static const uint64_t FIX2 = 0x3b8d12ac548a7c7aULL;
		static const uint64_t PCV1 = 0x3d84e1ac0dc82880ULL;
		static const uint64_t PCV2 = 0x0000000000000001ULL;
	};


	template<uint32_t MExp=19937>
	class SFMT
	{
	public:
		typedef SFMTParams<MExp> params;

	public:
		SFMT()
		{
			uint32_t seed = (uint32_t)std::time(nullptr);
			init(seed);
		}

		SFMT(uint32_t seed)
		{
			init(seed);
		}

		uint32_t gen_u32()
		{
		    if (m_idx >= params::N64)
		    {
		    	gen_rand_all();
		    	m_idx = 0;
		    }

		    uint64_t *psfmt64 = &m_status[0].u[0];
		    uint32_t r = psfmt64[m_idx++] & 0xffffffffU;
		    return r;
		}

		double gen_f64()
		{
		    if (m_idx >= params::N64)
		    {
		    	gen_rand_all();
		    	m_idx = 0;
		    }

		    double *dsfmt64 = &m_status[0].d[0];

		    union
		    {
		    	double d;
		    	uint64_t u;
		    } r;
		    r.d = dsfmt64[m_idx++];
		    r.u |= 1;

		    return r.d - 1.0;
		}

	private:
		void period_certification();

		void init(uint32_t seed);

		void gen_rand_all();

		void do_recursion(w128_t *r, w128_t *a, w128_t *b, w128_t *u)
		{
		    __m128i v, w, x, y, z;

		    x = a->si;
		    z = _mm_slli_epi64(x, params::SL1);
		    y = _mm_shuffle_epi32(u->si, 0x1b);
		    z = _mm_xor_si128(z, b->si);
		    y = _mm_xor_si128(y, z);

		    v = _mm_srli_epi64(y, params::SR);
		    w = _mm_and_si128(y, sse2_param_mask);
		    v = _mm_xor_si128(v, x);
		    v = _mm_xor_si128(v, w);
		    r->si = v;
		    u->si = y;
		}

	private:
	    __m128i sse2_param_mask;
	    __m128i sse2_int_one;
	    __m128d sse2_double_two;
	    __m128d sse2_double_m_one;

	    w128_t m_status[params::N + 1];
	    uint32_t m_idx;
	};


	// Typedefs

	typedef SFMT<19937> sfmt_19937;


	template<uint32_t MExp>
	void SFMT<MExp>::init(uint32_t seed)
	{
	    uint32_t *psfmt = &m_status[0].u32[0];
	    psfmt[0] = seed;

	    for (uint32_t i = 1; i < (params::N + 1) * 4; i++)
	    {
	        psfmt[i] = (uint32_t)1812433253UL
		    * (psfmt[i - 1] ^ (psfmt[i - 1] >> 30)) + i;
	    }

	    // initial mask

	    uint64_t *psfmt64 = &m_status[0].u[0];
	    for (uint32_t i = 0; i < params::N * 2; i++) {
	        psfmt64[i] = (psfmt64[i] & params::LOW_MASK) | params::HIGH_CONST;
	    }

	    period_certification();
	    m_idx = params::N64;

	    // setup const

	    sse2_param_mask = _mm_set_epi32(
	    		(int32_t)params::MSK32_3,
	    		(int32_t)params::MSK32_4,
	    		(int32_t)params::MSK32_1,
	    		(int32_t)params::MSK32_2);
	    sse2_int_one = _mm_set_epi32(0, 1, 0, 1);
	    sse2_double_two = _mm_set_pd(2.0, 2.0);
	    sse2_double_m_one = _mm_set_pd(-1.0, -1.0);
	}



	template<uint32_t MExp>
	void SFMT<MExp>::gen_rand_all()
	{
	    w128_t lung = m_status[params::N];
	    do_recursion(&m_status[0], &m_status[0], &m_status[params::POS1], &lung);

	    uint32_t i;
	    for (i = 1; i < params::N - params::POS1; i++)
	    {
	    	do_recursion(&m_status[i], &m_status[i], &m_status[i + params::POS1], &lung);
	    }

	    for (; i < params::N; i++)
	    {
	    	do_recursion(&m_status[i], &m_status[i], &m_status[i + params::POS1 - params::N], &lung);
	    }

	    m_status[params::N] = lung;
	}

	template<uint32_t MExp>
	void SFMT<MExp>::period_certification()
	{
		uint64_t pcv[2] = {params::PCV1, params::PCV2};
		uint64_t tmp[2];
		tmp[0] = (m_status[params::N].u[0] ^ params::FIX1);
		tmp[1] = (m_status[params::N].u[1] ^ params::FIX2);

		uint64_t inner = tmp[0] & pcv[0];
		inner ^= tmp[1] & pcv[1];
		for (int i = 32; i > 0; i >>= 1)
		{
			inner ^= inner >> i;
		}
		inner &= 1;

		// check OK
		if (inner == 1) return;

		// check NG
		if ((params::PCV2 & 1) == 1)
		{
			m_status[params::N].u[1] ^= 1;
		}
		else
		{
			for (int i = 1; i >= 0; i--)
			{
				uint64_t work = 1;
				for (int j = 0; j < 64; j++)
				{
					if ((work & pcv[i]) != 0)
					{
						m_status[params::N].u[i] ^= work;
						return;
					}
					work = work << 1;
				}
			}
		}
		return;
	}



}

#endif
